home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48hor1 / fft.src < prev    next >
Text File  |  1991-05-29  |  2KB  |  63 lines
  1. %%HP: T(3)F(,);
  2. @ by Jurjen NE Bos
  3. DIR
  4.   CST { { "FFT"
  5.     \<< 1 FFT
  6.     \>> } { "-FFT"
  7.     \<< -1 FFT
  8.     \>> } }
  9.   FFT   @ FFT routine.
  10.         @ Input:        2: Vector 1: 1 for FFT, -1 for inverse
  11.         @ Ouput:        1: transformed vector
  12.     \<< OVER SIZE 1 GET
  13.       IF OVER 0 <       @ Divide by length if inverse FFT
  14.       THEN ROT OVER / ROT ROT
  15.       END (0;6,28318530718) ROT * OVER / SWAP DUP
  16.       WHILE DUP 2 MOD NOT       @ Divide out factors of 2
  17.       REPEAT 2 /
  18.       END SWAP OVER / \-> \Gr r p       @ p: twopower, r: odd rest
  19.                                         @ rho: root of unity
  20.       \<< { r p } RDM
  21.         IF r 1 \=/      @ Compute regular r-point FFT p times simultaneously
  22.         THEN A\->L \Gr p * EXP 1 \-> m \Grp \Grn
  23.           \<< 1 r
  24.             START m 1 GET 1 2 r
  25.               FOR k \Grn * m k GET OVER * ROT + SWAP
  26.               NEXT DROP OBJ\-> DROP \Grp '\Grn' STO*
  27.             NEXT
  28.           \>> { r p } \->ARRY
  29.         END TRN CONJ    @ TRN does transpose & conjugate, so conjugate back
  30.         WHILE p 1 \=/   @ combine r-point FFTs to 2r-point FFTs
  31.         REPEAT OBJ\-> DROP 'p' 2 STO/ { p r } \->ARRY TRN A\->L
  32.         @ from here on, we work with conjugated numbers!
  33.           \Gr NEG p * EXP 1 \-> m \Grp \Grn
  34.           \<< { p r } \->ARRY TRN 1 r
  35.             FOR k m k GET \Grn * OBJ\-> DROP \Grp '\Grn' STO*
  36.             NEXT
  37.           \>> { r p } \->ARRY + LASTARG - \-> m
  38.           \<< OBJ\-> DROP m
  39.           \>> OBJ\-> DROP 'r' 2 STO* { r p } \->ARRY TRN
  40.           @ numbers are normal again
  41.         END { r } RDM
  42.       \>>
  43.     \>>
  44.   A\->L @ convert matrix to list of rows.  This doesn't use sigma+, because
  45.         @ this trick doesn't work with complex numbers :-(
  46.     \<< OBJ\-> OBJ\-> DROP { } \-> r c u
  47.       \<< 1 r
  48.         START c \->ARRY 'u' STO+
  49.         NEXT u
  50.       \>>
  51.     \>>
  52.   DFT   @ The regular Fourier Transform.  Extra short version.
  53.     \<< SWAP DUP SIZE 1 GET (0;6,28318530718) OVER / 4 PICK * \-> d x s q
  54.       \<< 0 s 1 -
  55.         FOR k '\GS(n=0;s-1;x(n+1)*EXP(q*k*n))' \->NUM
  56.         NEXT s \->ARRY
  57.         IF d -1 ==
  58.         THEN s /
  59.         END
  60.       \>>
  61.     \>>
  62. END
  63.